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An automatic numerical quadrature routine (ANQR) attempts to 

evaluate 

r b 

J a fMdX 

to absolute accuracy e, given only e, a, b, and a user-supplied subroutine 
which calculates f(x) for any x in [a,b]. An ANQR which guarantees 
success is impossible to construct, even disregarding the effects of finite 
computer precision, but the problem is nonetheless of interest. A reli- 
able and efficient ANQR is a necessary part of any mathematical sub- 
routine library. New single- and double-precision ANQRs, QUAD and 
DQUAD, have been constructed and tested. They are based on adaptive 
Romberg extrapolation, with cautious error estimation. An important 
practical feature is the automatic recognition of endpoint singularities, 
and a change of variable to handle them. QUAD and DQUAD also rec- 
ognize the presence of noise in the function being integrated, and limit 
the attempted accuracy accordingly. Since guaranteed ANQRs are 
impossible, extensive testing of DQUAD is presented to demonstrate 
its efficiency and robustness. Comparable testing is not available for 
competitive ANQRs, but performance on a standard set of test integrals 
is presented for DQUAD and nine other ANQRs. DQUAD is generally 
better. QUAD and DQUAD are written in PFORT, a subset of American 
National Standard (ANS) Fortran. Machine-dependent constants are 
obtained from the PORT library machine-constants programs. A por- 
table package of storage allocation routines is used. 

I. INTRODUCTION 

The development of automatic numerical quadrature routines 
(ANQRs) has been a popular research topic for many years (see refs. 1-6, 
8, 9, 11, 13, 14). An ANQR is a routine which attempts to calculate 

f(x)dx 
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with absolute error, or perhaps relative error, no larger than e, given e, 
a, b, and a procedure which calculates f(x) for any desired x in the in- 
terval [a,b]. It is assumed that no other information about the function 
/ is available. The problem is perhaps the more interesting for being an 
impossible one. Any numerical quadrature routine must estimate the 
integral by sampling the function / at a finite number of x's. A guaran- 
teed automatic integration algorithm is clearly impossible for general 
/, even for analytic /. For example, given any deterministic rule for nu- 
merical quadrature, one can readily find constants a and /3 so that the 
quadrature rule calculates 



^So e " 



to be close to zero. (Choose a to be large and positive, and /3 to be between 
sampling points.) 

Although the general problem is impossible, one feels that an ANQR 
which works for "reasonable" functions should be feasible, and much 
work has been directed at this goal. There has been great confusion and 
difficulty in comparing the various candidates for ANQRs, partly because 
the domain of the problem is undefined; a reasonable definition of a 
"reasonable" function is itself difficult. 

In constructing an ANQR, an author is forced to make decisions about 
the class of "reasonable" functions, in effect to define what is a "rea- 
sonable" function. These decisions strongly affect the efficiency and 
robustness of the ANQR. For example, to avoid completely missing an 
isolated peak in f(x), the interval [a,b] must be sampled finely. However, 
a fine sampling is inefficient for easy functions. Another example is a 
function which is flat over 99 percent of the interval, and which has 200 
oscillations in the remaining 1 percent. If an ANQR is able to distinguish 
this function from one which is merely noisy over 1 percent of the in- 
terval, the ANQR is likely to be inefficient on easy functions and very 
inefficient on noisy functions. An ANQR which gives up relatively quickly 
on this function, calls it noisy, and returns an error message, may be 
preferable, especially since many such functions are the result of a user's 
programing errors. 

A compromise strategy, used by QUAD, is to isolate all assumptions 
about the "reasonable" class of functions in a few parameters. Default 
values of these parameters can be chosen which will be suitable for most 
users. More knowledgeable users can use other values. With the default 
values, QUAD strikes what the author considers to be the proper balance 
between efficiency and robustness. 

Since no a priori information about f(x) is available, 

f(x)dx 
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must be evaluated by sampling f in [a, b]; the error in the calculated in- 
tegral is usually estimated by comparing two or more calculated values 
for the integral. ANQRs typically have a sequence of quadrature rules 
Q n , depending on a, b, and the function /, such that 



J. b 
f(x) dx 
a 



if the calculations are done in infinite precision, and if / is at least 
piecewise continuous. Most ANQRs have no better error estimation 
procedure than to accept Q n whenever \Q n -i ~ Qn\ < «, a procedure 
fraught with danger. QUAD has a much more stringent error estimation 
procedure, described in Section II. 

Many functions to be integrated are easy to integrate over some parts 
of the interval and difficult over other parts. It is frequently more effi- 
cient to sample more densely in the difficult regions, if possible. ANQRs 
which attempt to do this are called adaptive — the points at which / is 
sampled depend on the function being sampled. An adaptive ANQR must 
include some strategy for how to concentrate the sampling points. Es- 
sentially all competitive ANQRs are adaptive. 

The usual adaptive procedure is to integrate an interval with quad- 
rature rules Q n for n » 1, 2, . . . , N, where N is fixed. Q n may be, for ex- 
ample, Simpson's rule with 2 N intervals, or Gauss-Legendre quadrature 
with n sampling points. If convergence has not been obtained, the in- 
terval is divided in half, and each half considered separately. For effi- 
ciency, one wants quadrature rules for which all sampling points for the 
whole interval are also used for the half- intervals. If the value of N used 
depends on the results Q n for n < N, the method is sometimes called 
doubly -adaptive. 

Most ANQRs do not do well on integrals with endpoint singularities, 
but users' integrals are frequently of this type. QUAD has a provision 
for recognizing endpoint singularities and for making a change of variable 
to facilitate the integration. This feature also works well on another 
important class of functions, those decaying steeply away from one or 
both ends of the integration interval. This automatic change of variable 
technique is a significant improvement over previous ANQRs. 

Most ANQRs cannot cope with noisy functions; if there is too much 
noise in /, most ANQRs fail in an unpleasant, uneconomical way. Con- 
vergence will be at best very slow, so that the ANQRs will stop only when 
their predefined limit on calls to the function evaluation procedure has 
been exceeded, with no indication that the problem is noise rather than 
a noise-free but unruly function /. QUAD recognizes noisy functions, sets 
a warning flag, and integrates only to an accuracy commensurate with 
the estimated noise. 

Finally, there is a large difference between an algorithm for numerical 
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quadrature and a properly-written ANQR suitable for a program library. 
Provision must be made, for example, to stop trying to integrate a 
function if it has been sampled more than some user-defined number 
of times. The finite machine precision of the computer involved must 
be taken into account. Temporary storage must not be allowed to over- 
flow. Provision for error returns must be made. 

The basic idea behind QUAD is adaptive Romberg extrapolation 5 
combined with cautious error estimation 9 . The first such combination 
was the program cadre, written by deBoor 6 . cadre and QUAD are su- 
perficially similar, but differ in almost every detail. The major im- 
provements incorporated in QUAD include the following, which will be 
covered fully in Section II. 

(i) Noise. QUAD detects noisy functions and quits gracefully. 

(ii) Endpoint singularities. QUAD detects singularities in f(x) at the 
endpoints, a and b, and automatically makes a change of variable to 
reduce the strength of the singularity. 

(Hi) Mesh sequence. QUAD uses the mesh sequence 1, 2, 3, 4, 6, 8, 12, 
16, ... , instead of 1, 2, 4, 8, 16, ... , giving a higher effective order of 
convergence. 

(iv ) Portability. QUAD is written in PFORT, 15 a portable subset of ANS 
Fortran. Machine-dependent quantities are defined with the PORT 7 
machine constants. A portable Fortran stack 7 is used for temporary 
storage. 

Section II discusses the algorithm of QUAD and DQUAD more fully. 
Section III compares the performance of DQUAD and nine competitive 
ANQRs on a standard set of test integrals, and also presents the results 
of some more serious testing of QUAD. Section IV discusses the imple- 
mentation of QUAD, including portability considerations. 

II. QUAD 

2.1 Romberg extrapolation 

QUAD is based on Romberg extrapolation of the composite trapezoidal 
rule. 5 The formulas are standard, but will be repeated here for com- 
pleteness. Let Mi, M2, ... be an increasing sequence of positive integers, 
and let hi = (b — a) Mi. Then the composite trapezoidal approximation 
to 



= J> 



b 

x) dx 



IS 



T(hi) = 1 -h i \f{a) + f{b)] + hi n Zf(a + m/i f ) 

2, m=l 
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If/ has 2fc + 1 continuous derivatives in [a,b], then the Euler-Maclaurin 
sum formula shows that 

T(hi) = I+ E c m hf m + 0(/iP +1 ) 

m = l 

where the c m depend only on a, b, and /, not on hi. A higher-order, al- 
though not necessarily more accurate, estimate may be obtained by 
combining two trapezoidal estimates via Richardson extrapolation, 3 
eliminating the C\hf term. Let Tft ] = T(hi). 

T (l) _ rp(2) ■ T o ~ ^6 

= 1 + 0{h\hl) 

Still higher-order estimates may be generated recursively. The general 
formula for generating T^ is 

*r(i+l) _ T(') 
rp(i) _ r(i'+l) i 1 «-l L K-l 

n ~ n "' + hf/hh k - 1 

It is customary to think of the T-values as a table, viz. 
y(i) 

7f > T[ l) 

T h) T j2) T jl) 

Ti 4) T} 3) rp T^ x) 

T ^5) T J4) T ^3) T £2) T U) 



The classical Romberg method uses the sequence 1, 2, 4, 8, ... for the 
Mt's. 

Since at any time the right-most element, the "tip," of the T table is 
of highest order, it is expected to be most accurate, and frequently is so. 
Early Romberg programs 1 tested for convergence only by checking 
successive tip elements of the table. There are several arguments against 
this practice. Firstly, for the tip there is no way to obtain an error esti- 
mate with any theoretical foundation. Secondly, highest-order is not the 
same as most accurate. Even for analytic functions, if the step size used 
at the beginning of the T table is too large, the tip of the table may not 
be the most accurate value. Thirdly, for functions f(x) which do not have 
enough derivatives, the tip of the T table is not of higher order than the 
elements to the left. For example, if f(x) = x a , for < a < 1, the low- 
est-order term in each column of the feth row is 0{h a+1 ). In practice it 
is frequently found that lower-order columns are more accurate than 
the tip element. 
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2.2 Cautious error estimation 

The idea of cautious error estimation comes from Lynch. 9 It is a simple 
and seemingly unobjectionable idea, but is not adopted by most authors 
of ANQRs. A cautious error estimation procedure believes an error esti- 
mate only if there is some evidence that the convergence rate of suc- 
cessive quadrature rules is close to the theoretical rate. Cautious error 
estimation is particularly easy for Romberg extrapolation. The two were 
first combined by deBoor. 6 quad's version of catious error estimation 
is similar in spirit to deBoor's, but is more cautious and is different in 
all details. 

It has been proven 1 that, if/ is merely Riemann-integrable, each col- 
umn of the T table converges, as does each diagonal. If / has enough 
continuous derivatives, 

n»=i+o(hf...hf +k ) 

These theoretical results provide a basis for cautious error estimation. 6,9 
Lynch's suggestion was to consider three successive trapezoidal rule 
estimates, and form the ratios 

„., _ n> - n*° 

TH +1) — Tl' +2) 

hf-h? +1 + 0[(hf + hl +1 )*] 

h? +1 -hf +2 + 0[(hf +1 + hf +2 )*] 

If the step sizes are small enough, the higher-order terms are small 
compared to the second-order terms, and 



ftj« 



hf-hhi 
hi+i — hf+z 



The calculated R^ l) being close to this theoretical value is good evidence 
that the convergence rate of the column is proper, and that the error in 
Tb l) is dominated by C\hf. Then the Runge estimate of the error, 3 

Tff+2) _ *T(j+l) 



hf +l /hf +2 -l 

is likely to be a good estimate. If the calculated R^ l) is not close to the 
theoretical value, the Runge estimate of the error is likely to be an un- 
derestimate. 

Similar calculations are done in higher columns. The general formula 
forflifMs 

R p m n ] -n +i) 



jffl) _ T j.+2) 
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As all the /i's approach zero, 

^ ~ i 9 i 2 i' 

"i+fc+1 "i+l — n T+k+2 

The Runge estimate of the error in Tjt? +2) is 

|Tjf + »-/|* 7^7 if-: 

hi 2 /hi+k+i 2 - 1 

A more conservative estimate of the error in Tj/ +2) , 



|7f+2)_ 7 | «|7t+ 2) -7t' +1) | 

is often used. 

QUAD calls column k asymptotic if Rtf is close enough to the theo- 
retical value; the tolerance is 5 percent of the theoretical value for k = 
0, 10 percent for k ■ 1, 15 percent for k = 2, and so on. Column fc is al- 
most asymptotic if Rf is between 0.25 and 4.0 times the theoretical 
value, except for column 0, where the criteria are 0.75 and 1.25. 

If columns through k are asymptotic, QUAD believes the Runge es- 
timate for the error in the kih column. If columns through k - 1 are 
asymptotic and column k is almost asymptotic, QUAD believes the 
conservative estimate for the error in the fcth column, but no higher 
columns are believed. (The only exception is that, if the column is only 
almost asymptotic, the next column is believed if it itself is asymptot- 
ic.) 

This describes the basic cautious error estimation procedure for QUAD. 
There are a few more details, however. QUAD does not believe any answer 
based on less than two extrapolations, or five sampling points, per in- 
terval. If two successive entries in a column give the same value to within 
a few rounding errors, as occurs when integrating a constant function 
or in doing very accurate integration, then the column does not appear 
to be asymptotic. The conservative error estimate for the column is be- 
lieved anyway. If an interval has a singularity, either real or due to 
rounding errors or truncation errors in the function subprogram for f{x), 
no column will appear asymptotic. Then a nonasymptotic answer in the 
first column will be accepted after several extrapolations, with the very 
conservative error estimate 

\Tg+* ~I\ * 2|n ,+2) - Tjj +l) \ + 2|T# +1) - r °l 

Finally, at any stage one column of the T table has only two entries in 
it, and cannot be judged to be asymptotic or nonasymptotic. The con- 
servative error estimate is accepted for such a column if the previous 
column is asymptotic. 
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2.3 Step size sequence 

The above discussion applies to any sequence of step sizes. The clas- 
sical Romberg sequence uses step sizes (6 — a)/n, with n ■ 1, 2, 4, 8, 16, 
. . . , halving the step size and doubling the number of sampling points 
at each new extrapolation. Several alternative sequences have been 
suggested which do not cause the number of sampling points to rise so 
rapidly. QUAD uses n's of 1, 2, 3, 4, 6, 8, 12, . . .; another reasonable pos- 
sibility is 1, 2, 3, 4, 5, 6, 8, 10, 12, These sequences double the number 

of sampling points every second and third extrapolation, respectively. 
The classical sequence has the advantages that the bookkeeping is very 
easy and that all old sampling points are reused if the interval is divided 
in half. The latter is essential for efficiency. 

QUAD'S sequence uses fewer sampling points to get the same accuracy, 
as suggested by Bulirsch and Stoer. 3 The bookkeeping is more compli- 
cated than for the classical sequence. For example, it is only convenient 
to divide the interval in half after the fourth, sixth, eight, . . . extrapo- 
lations if all the old sampling points are to be reused. 

2.4 Adaptive procedure 

For an adaptive Romberg extrapolation routine, it is necessary to 
decide when to do another trapezoidal rule and another extrapolation, 
and when to divide the interval. The minimum number of extrapolations 
for QUAD is 4, using step sizes (6 — a) down through (b — a)/6, and a 
total of 9 sampling points. This default lower limit may be raised by the 
user (see Section IV). Because of roundoff, unlimited extrapolations are 
impractical — the highest-order columns will not be asymptotic and will 
not be believed. The maximum number of extrapolations allowed by 
QUAD is 6; DQUAD allows 8. The default limit may be changed by the user 
(see Section IV). 

If the requested error tolerance for an interval has not been achieved 
after 4 extrapolations, QUAD goes on to 5 and 6 extrapolations if the first 
column is asymptotic; after 6, it goes on if the second column is also as- 
ymptotic. This procedure is biased in favor of doing more extrapolations, 
and trying to get higher-order convergence, for smooth functions. 
Functions which are not smooth, or which do not appear asymptotic 
because of too large a step size, have the interval divided instead of 
having more extrapolations done. If QUAD decides to divide an interval, 
the lower half is stacked, and the upper half is attempted next. 

2.5 Change of variable 

Functions with singularities are expensive to integrate without special 
methods. Since the interval containing the singularity will have a 
nonasymptotic T table, convergence will be limited to the first column. 
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For example, it can be shown 10 that for 

/= C 1 f(x)dx= Cx«g(x)dx (1) 

Jo Jo 

where g is smooth, and /(0) is set equal to zero, 

Ti(h) = 1+ t c m h 2m + £ d m h«> + « (2) 

The dominant error term for the first column of the T table is likely to 
be the d\h l+a term, so convergence is slow. If /(0) is not zero, another 
infinite sum is added to (2), like the second sum, but with a = 0. 

For such an endpoint singularity, the error is of a simple form. It is 
feasible to recognize this type of singularity in the same way that the 
cautious error estimation procedure recognizes asymptotic, or h 2m be- 
havior. De Boor 6 does exactly this, estimates an a, and then extrapolates 
using eq. (2). The success of this procedure depends critically on how 
accurately a can be estimated. If /(0) is not set equal to zero, de Boor's 
method will not work well. For logarithmic singularities, the error ex- 
pansion corresponding to eq. (2) is more complicated; de Boor makes 
no attempt to recognize logarithmic singularities. 

After recognizing an endpoint singularity, QUAD uses a different 
procedure. Suppose that the integral is as above, where g is well-behaved. 
Then the leading error term is 0{h 1+a ) if — 1< a < 1. In the second and 
higher columns, the 0(h 2 ) term is gone, so the leading term is 0(h 1+a ) 
for —1 < a < 3. QUAD looks at ratios of T table entries in the second and 
third columns to recognize x" behavior, and estimates the value of a. 
QUAD then makes a change of variable x = u n , where n is the closest 
integer to 6/(1 + a), giving for eq. (1) 

f nu n ~ 1 f(u n ) du = f nu n ^ +a) - x g{u n ) du 

(The change of variable is somewhat more complicated if the limits of 
integration are not and 1.) The new integral has a singularity of the 
form u$, where /3 is between 4.5 and 5.5, if a is close to the true a. The 
singularity in the transformed integral is lessened, allowing convergence 
in the second or third columns of the T table. Convergence is likely to 
be much quicker. For rapid convergence, the method does not rely on 
the estimated a being close to the true a or upon eq. (2) holding, or indeed 
on there being any singularity at all at the endpoint. 
Steeply decaying integrands such as 

100 Ce-™ x dx 

look like step functions when coarsely sampled. A step function is an x° 
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singularity, since /(O) = 0, so a change of variable is made with n = 6. 
Singularities of the form x& log (x) will look sufficiently like x a singu- 
larities for some a, so that the transformation will be made. (It is not 
obvious that this is true, but tests have strongly indicated that it is.) 

If a is close to — 1, n can become large. QUAD requires n to be less than 
a maximum value determined by the precision of the computer; this 
value was 22 for tests reported in Section III. The change of variable is 
not made if a is less than -0.99. 

To facilitate the change of variable, quad starts by dividing the in- 
terval [a,b] into three equal intervals, and reverses the upper third. 
(Three is the default number, and may be changed by the user — see 
Section IV.) On the lower third and the reversed upper third, a left-hand 
endpoint singularity is recognized by a pattern of "fail on whole interval, 
succeed on right half- interval" twice in a row, and is followed by the 
estimation of a. If the two estimated values for a from the second and 
third columns of the T table do not agree to within 0.1, no change of 
variable is made. No change of variable is attempted except at the two 
endpoints of the originial interval. 



2.6 Noisy functions 

All procedures for evaluating f(x) are inherently noisy, since they are 
implemented on finite-precision machines. The value returned is not 
the exact f(x), but/(x) [1 + r\(x)] + r2(x), where r\(x) and r2(x) are noise 
functions. Ideally, r\ could always be no larger than a few rounding errors, 
and r 2 could be no larger than a few times the smallest positive ma- 
chine-representable number. Noise of this size should not affect the 
performance on an ANQR unless e is very small, of the order of r 2 (x) or 
f(x) ri (x). 

Protecting against this magnitude of noise is quite easy, although few 
ANQRs bother to do so. QUAD estimates a priori the sizes of r\ and r 2 , 
based on the machine precision, and requires all error tolerances to be 
at least as large as the estimated rounding error. 

Protecting against significantly larger noise is more difficult. A suc- 
cessful ANQR should recognize the presence of noise, estimate its mag- 
nitude, and evaluate the integral in a "reasonable" number of function 
evaluations with an accuracy which is "nearly" as good as possible. 
(There is of course a trade-off between "reasonable" and "nearly.") If 
typical values of r%{x) or of f{x)r\(x) are much larger than e/\ b — a |, most 
ANQRs will fail in an unpleasant, uneconomical way. Convergence will 
at best be very slow, so that the ANQRs will stop only when their prede- 
fined limit on calls to the function evaluation procedure has been ex- 
ceeded, with no indication that the problem is noise rather than a 
noise-free but unruly function /. 
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SQUANK 11 makes a reasonable effort at recognizing noise. However, 
it attributes any nonstandard behavior to noise, so that some unruly but 
noise-free functions are called noisy. 

Qualitatively, one may say that a function is noisy if, on a "sufficiently 
small" interval, the values of the samples of/ are "not smooth enough." 
An algorithm consists of the defining of "sufficiently small" and "not 
smooth enough," followed by estimation of the magnitude of the noise 
and by further action to avoid using an excessive number of function 
values. 

In QUAD, no answer is believed unless the function has been sampled 
with adjacent samples no farther apart than h s \ b - a | ; h s is supposedly 
small enough so that all structure may be seen by sampling with this 
spacing (h s is parameter HSAMPL of Section IV). The default value of 
h s is Vg. Noise is not estimated unless adjacent samples are no farther 
apart than h n \b - a | , with h n = h s /S2. Choosing h n smaller would re- 
quire more function values; choosing h n larger would increase the risk 
of calling a function noisy when it is noise-free but rapidly varying. 

When a function is noisy, QUAD will usually fail on large intervals, and 
then attempt smaller and smaller subintervals. When integration on a 
subinterval has failed, and adjacent samples are no farther apart than 
h n \b — a\, the noise in / is estimated. First, the second differences of the 
samples of/ on the subinterval are formed, e.g. f{x) — 2f(x + h) + f(x 
+ 2h). If the sequence of second differences has no more than two sign 
changes, noise is not assumed to be present. If there are three or more 
sign changes, noise is assumed to be present, since the function has too 
much structure over too small an interval, and the estimated answer and 
error for that subinterval are accepted as being as good as possible. 

When noise has thus been found to be present, the second differences 
are assumed to be essentially all noise, and the magnitude of the noise 
is estimated as the average of the absolute values of the second differ- 
ences. All succeeding subintervals are attempted with accuracy not ex- 
ceeding the estimated noise magnitude times the length of the subin- 
terval. If other subintervals are found to be noisy, the largest noise 
magnitude is used. 

2.7 Error allocation 

QUAD attempts to integrate the upper one-third of [a,b] with error 
tolerance t/3. Then the following procedure is used to assign an at- 
tempted accuracy for each interval. When integration on an interval is 
attempted with error tolerance e and fails, the upper half is attempted 
with error tolerance «o/2. When integration on an interval succeeds, the 
absolute value of the estimated error is added to a running sum, and the 
top interval is the stack is attempted. If the top interval is of length bx, 
the total length of intervals remaining to be integrated is Xi, and the sum 
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of the absolute values of the error estimates so far is €1, then the error 
tolerance assigned to the top interval in the stack is (e — ei) bx/x\. 
However, the error assigned to any interval is required to be at least as 
large as e/1000. 

If an answer is returned for an interval with an error estimate which 
is larger than the requested error, but less than the estimated roundoff 
or noise, that answer and error are accepted, and a warning flag is set. 



III. TESTING AND COMPARISON OF ROUTINES 

Testing is necessary to evaluate the efficiency and robustness of an 
ANQR. Typically the proposer of an ANQR generates an algorithm, codes 
a simple program, and tests the routine on a few integrals of the pro- 
poser's own choosing. It is not unusual for all the test integrals to be done 
well by the ANQR. As a reqult, the prospective user has no way of eval- 
uating the quality of the ANQR without performing extensive testing. 

Some improvement was evidenced in the work of Kahaner, 8 who tested 
many ANQRs on the same set of 21 test integrals. The same set was used 
by de Boor 6 for testing his ANQR. At least three of the 21 are not appro- 
priate test integrals, however, because the results are algorithm-de- 
pendent in an unrepresentative way. 

Two examples will make this clear. First consider 



£ 



i 

f(x) cos (ccirx) dx 





where f(x) is any smooth function. QUAD, which divides according to the 
1> %« %i Vii % sequence, will fail for a near 36, but not for a near 32. [For 
a near 36, the regular sampling procedure of QUAD samples only near 
the peaks of the cosine, so the integrand looks like f(x).] CADRE, 6 which 
divides according to the 1, V2» l U> Vs sequence, will fail for a near 32, but 
not for a near 36. A single test integral with a large a may not compare 
ANQRs fairly. Test integrals 13 and 17 of Kahaner (see Table I) are of 
this type, each having about 50 full cycles. 
Second, consider 



-*$>- 



for a large and positive. Depending on the choice of ft this can be either 
easy or hard for a particular ANQR. If the peak comes sufficiently near 
a sampling point, adaptive ANQRs can zero in on the peak and integrate 
it accurately, although many sampling points will be necessary. If the 
peak does not come sufficiently near a sampling point, the integrand 
looks like zero, and the ANQR fails. For proper comparison of ANQRs, any 
single a is insufficient. Test integral 21 of Kahaner is of this type. 
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Table I — Kahaner's 21 Test Integrals 



No. 


a 


b 


answer 


/(x) 


Easy 
12 





1 


+0.7775046341 


x/(e* - 1) 


1) 





1 


+0.3798854930 


1/(1 + e x ) 


1 





1 


+ 1.7182818284 


e x 


10 





1 


+0.6931471806 


1/(1 + x) 


4 


-1 


1 


+0.4794282267 


0.92 cosh (x) -cos(x) 


8 





1 


+0.8669729873 


1/(1 + x 4 ) 


5 


-1 


1 


+ 1.5822329637 


l/(x 4 + x 2 + 0.9) 


20 


-1 


1 


+ 1.5643964441 


1/U 2 + 1.005) 


Steeply 










decaying 










15 


n 


10 


+1.0000000000 


25c- 251 


14 





10 


+0.5000002112 


-v/50exp(-507rx 2 ) 


16 





10 


+0.4993638029 


50/[p(l + 2500x 2 )] 


Singular 










6 





1 


+0.4000000000 


x 3/2 


3 





1 


+0.6666666667 


x 1/2 


2 





1 


+0.7000000000 


0, x < 0.3; 1, x > 0.3 


7 





1 


+2.0000000000 


0,x = 0;x" 1/2 ,x >0 


19 





1 


-1.0000000000 


0,x = 0;ln(x),x>0 


Oscillatory 
18 





p 


+0.8386763234 


cos [cos x + 3 sin x + 2 cos 2x + 
3 sin 2x + 3 cos 3x] 


9 





1 


+0.4794282267 


2/[2 + sin (10px)] 
50[sin(50px)/(50px)] 2 


17 


0.01 


1 


+0.1121395696 


13 


0.1 


1 


+0.0090986453 


sin (100px)/px 


Isolated peak 
21 





1 


+0.2108027354 


sech 2 [10(x -0.2)] + sech 4 [100(x - 
0.4)] + sech 6 [1000(x - 0.6)] 


Note: p = 3.14 


59, P- 


= 3.1415927. 





Kahaner did not test noisy functions, and did not ask for impossibly 
small error tolerances. 



3. 1 Testing on the Kahaner 2 1 

The Kahaner 21 are listed in Table I. They have been grouped ac- 
cording to type; within groups they are ordered approximately by dif- 
ficulty. Tables II, III, and IV summarize the results often ANQRs on the 
21 test integrals, for requested error tolerances 10 -3 , 10 -6 , and 10 -9 . In 
each table, the first column is the integral number; succeeding columns 
are the number of sampling points used by each ANQR. An F indicates 
that the ANQR was unsuccessful, and an asterisk that the ANQR used 
the fewest sampling points of any successful routine. (For these tables 
only, "successful" means that the true error of an integration is no more 
than 20% higher than the requested error, since some ANQRs failed by 
a small amount.) Columns labeled ROMB through RBUN are based on 
the number of sampling points reported by Kahaner 8 for seven of his 
highest-quality ANORs, and were obtained on a CDC 6600. Column 
CADRE is from de Boor, 6 and the integrals were performed on a CDC 6500. 
Columns QSUBA and DQUAD were computed especially for this com- 
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Table II — Number of sampling points used by each anqr; 
attempted absolute accuracy 10 -3 



Inte- 






















gral 


ROMB 


SIMPSN 


SQUANK QNC7 QNClO QABS 


RBUN 


CADRE 


QSUBA DQUAD 


12 


17 


19 


9 


25 


37 


13 


5* 


9 


7 


13 


11 


17 


19 


9 


25 


37 


13 


5* 


5* 


7 


13 


1 


17 


19 


9 


25 


37 


13 


5* 


9 


7 


13 


10 


17 


19 


9 


25 


37 


13 


5* 


9 


7 


13 


4 


17 


19 


9 


25 


37 


13 


5* 


17 


7 


13 


8 


17 


19 


9 


25 


37 


13 


5* 


9 


7 


15 


5 


17 


19 


9* 


25 


37 


13 


11 


33 


15 


17 


20 


17 


31 


9* 


25 


37 


13 


11 


17 


15 


19 


15 


513 


103 


53 


85 


109 


85 


527 


88 


63 


52* 


14 


1025 


103 


49* 


97 


127 


85 


51 


62 


3F 


52 


16 


2049 


115 


53 


121 


163 


109 


87 


81 


127 


48* 


6 


17 


19 


9 


25 


37 


13 


5* 


9 


7 


15 


8 


65 


55 


9F 


49 


55 


77 


211 


17 


15* 


17 


2 


257F 


115 


29F 


121 


163 


141 


271 


53* 


771 


85 


7 


8193F 


235 


105F 


241 


361 


133F 


211 


33 


517 


26* 


19 


4097 


175 


45 


217 


307 


181 


211 


137 


31 


28* 


18 


129 


139 


53* 


85 


73 


77 


39F 


107 


63 


61 


9 


33* 


163 


81 


97 


145 


149 


79 


183 


127 


101 


17 


1025 


151 


57* 


165 


307 


149 


109 


512 


255 


185 


L3 


1025 


19F 


429 


49F 


865 


573 


533 


1028 


255* 


381 


21 


4097 


127 


17F 


97 


127 


77 


65 


108 


333F 


49* 



* = best of all successful results 
F = failure 



parison, on a Honeywell 6070. Double precision was used so that the 
relative machine precision would be comparable to that of the CDC 
machines. QSUBA 14 has provision only for relative error; for these tests 
a relative error was requested which gave the appropriate absolute error 
request. 

It is important to notice that some of the failures are due to an ANQR 
deciding to stop because of excessive sampling; these failures are far less 
reprehensible than the others because an error return could have been 
made, and an incorrect answer rejected. For RBUN through ROMB, this 
information can only be inferred since Kahaner did not list any error 
returns. No such failure occurred for CADRE or DQUAD; QSUBA has no 
built-in maximum. 

ROMB 1 is a standard Romberg extrapolation routine, using the stan- 
dard 1, V2, V4, Vs sequence. It is not adaptive. It stops, apparently, after 
8193 sampling points of f(x). Thus only one of its failures is serious. 
ROMB requires at least 17 points before believing any answer. 

SIMPSN 5 and SQUANK 11 are adaptive Simpson's rule routines, ap- 
parently with cutoffs at 5003 and 5001 points, respectively. These rou- 
tines are decent at low accuracies, but are not of high-enough order to 
be competitive at high accuracy. SQUANK also assumes that an improper 
convergence rate is due to noise in the function, rather than to a singu- 
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Table III — Number of sampling points used by each anqr; 
attempted absolute accuracy 10 -6 



Inte- 






















gral 


ROMB 


SIMPSN SQUANK 


QNC7 


QNC10 


QABS 


RBUN 


CADRE 


QSUBA DQUAD 


12 


17 


19 


9 


25 


37 


13 


5* 


9 


7 


13 


11 


17 


19 


9 


25 


37 


13 


11 


9 


7* 


13 


1 


17 


55 


17 


25 


37 


13 


21 


17 


7* 


19 


10 


17 


55 


21 


25 


37 


13* 


31 


17 


15 


19 


4 


17 


55 


25 


25 


37 


25 


5F 


33 


15* 


19 


8 


33 


67 


29 


25 


37 


25 


41 


17 


15* 


21 


5 


65 


163 


65 


49 


37 


49 


59 


49 


31* 


41 


20 


65 


163 


49 


49 


37 


49 


49 


33 


31* 


33 


15 


2049 


343 


213 


133 


145 


133 


4117 


140 


63* 


98 


14 


2049 


331 


169 


133 


163 


109 


91 


89 


3F 


86* 


16 


8193F 


511 


273 


181 


181 


145 


141 


145 


255 


96* 


6 


129 


91 


29* 


61 


73 


65 


383 


65 


31 


34 


3 


4097 


199 


105 


157 


217 


145 


423 


33 


63 


32* 


2 


8193F 


235 


29F 


241 


361 


261 


271 


119* 


3351 


125 


7 


8193F 


1027F 


1153F 


241F 


361F 


89F 


587F 


129 


5925 


40* 


19 


8193F 


499 


257 


241 


361 


105F 


403 


233 


795 


38* 


18 


257 


547 


301 


181 


199 


205 


195 


177 


63* 


117 


9 


129* 


871 


377 


289 


397 


313 


267F 


409 


255 


285 


17 


2049 


2275 


697F 


385F 


1009 


829 


697 


1237 


255* 


547 


13 


2049 


19F 


2549 


1525 


1639 


1449 


2383 


1449 


255* 


757 


21 


8193F 


691 


185F 


205F 


253F 


197F 


327* 


189F 


525F 


127F 



* = best of all successful results 
F = failure 



larity, and so quits early on some of the test integrals. SIMPSN requires 
at least 19 points, and squank 9. 

QNC7 and QNClO 8 are adaptive Newton-Cotes routines, with 7- and 
10-point rules, respectively. (QNClO was called QUAD in Ref. 8.) They 
performed quite well, failing only on some of the most difficult integrals. 
QNC7 requires at least 25 points, and QNClO at least 37, somewhat ex- 
cessive for the easiest integrals. 

QABS 13 combines Romberg and Curtis-Clenshaw quadrature, and 
performed quite well, failing only on some of the most difficult integrals. 
It requires at least 13 points. 

RBUN 4 is an adaptive Romberg extrapolation routine, using the 
standard sequence. It apparently has a cutoff at 5001 points. RBUN re- 
quires at least 5 points, and seems to be somewhat unreliable. 

Kahaner recommended any of QNC7, QNClO, or QABS as a library 
routine. He did not have CADRE, QSUBA, or DQUAD available to test. 

CADRE 6 is more recent than the ANQRs just discussed. It uses a version 
of cautious Romberg extrapolation based on the standard sequence. It 
also includes provision for recognizing a singularity of the form x a , es- 
timates a numerically, and extrapolates using the estimated a. CADRE 
requires at least 5 points. On the test integrals, it seems somewhat more 
efficient than QNC7, QNClO, and QABS on nonsingular integrals and 
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Table IV — Number of sampling points used by each anqr; 
attempted absolute accuracy 10~ 9 



Inte- 






















gral 


ROMB 


SIMPSN SQUANK QNC7 


QNC10 


QABS 


RBUN 


CADRE 


QSUBA 


DQUAD 


12 


17 


55 


33 


25 


37 


13* 


39 


17 


15 


19 


11 


33 


151 


33 


25 


37 


25 


39 


17 


15* 


19 


1 


17 


163 


65 


25 


37 


25 


73 


17 


15* 


23 


10 


65 


271 


97 


37 


37 


49 


123 


33 


15* 


25 


4 


33 


331 


105 


25 


37 


85 


139 


33 


15* 


25 


8 


65 


463 


149 


73 


73 


97 


129 


65 


31* 


45 


5 


129 


487 


289 


97 


73 


181 


239 


129 


31* 


73 


20 


129 


487 


249 


97 


73 


145 


185 


129 


31* 


49 


15 


4097 


1483 


1145 


241 


217 


281 


5001F 


215 


127* 


138 


14 


4097 


1123 


797 


241 


253 


245 


259 


202 


3F 


154* 


16 


8193F 


2467 


1649 


397 


343 


397 


1435 


337 


255 


192* 


6 


2049 


427 


161 


133 


163 


137 


1423 


529 


63 


46* 


a 


8193F 


883 


513 


289 


361 


289 


1595 


129 


255 


42* 


2 


8193F 


235 


29F 


241F 


361F 


381 


271 


173 


5931 


165* 


7 


8193F 


4279F 


5001F 


589F 


685F 


89F 


2467F 


625 


11325 


70* 


19 


8193F 


2203 


1969F 


421F 


415F 


89F 


1571 


369F 


3495 


80* 


18 


513 


2923 


1589 


409 


343 


589 


753 


417 


127* 


217 


9 


257* 


3967 


2525 


697 


757 


893 


883 


785 


765 


473 


17 


4097 


5003F 


5001F 


1345F 


1999 


2025 


2741 


2329 


255* 


1109 


13 


4097 


5003F 


5001F 


3073 


2773 


3197 


5001F 


3505 


255* 


1161 


21 


8193F 


3751 


1657 


709 


685 


633* 


1079 


661 


827F 


261F 



* = best of all successful results 
F = failure 



much more efficient on singular ones. In addition, the cautious extrap- 
olation means that cadre's error estimation procedure has some ra- 
tionale behind it, and CADRE is more robust than the aforementioned 
routines. However, CADRE is difficult to understand and to maintain, 
since its style is the antithesis of structured programming. It is one large 
program, with no subprograms, but with a liberal and unstructured use 
of GOTOs. 

QSUBA 14 uses a series of 8 whole-interval quadrature rules of increasing 
order, starting with the 1- and 3-point Gauss-Legendre rules. Succeeding 
rules are constructed to be of as high order as possible, consistent with 
using all the previous sampling points. The highest-order rule uses 255 
points and is of order 383. If convergence is not obtained after 8 rules, 
the interval is divided in half, and each half considered anew. Unlike all 
the other ANQRs under consideration, all function values must be dis- 
carded, since none of the sampling points on the half intervals coincides 
with any on the full interval. QSUBA works well on any integral which 
can be integrated without dividing the interval, and poorly on integrals 
which require dividing. It is especially good on easy and oscillatory in- 
tegrals, since a high-order rule is generally used. QSUBA uses at least 3 
points, which is somewhat unsafe, but has no maximum built in — it goes 
on forever, if necessary. 
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DQUAD takes at least 13 points. It fails only on number 21, for high 
accuracy, missing the narrow peak at x = 0.6. The integral which is the 
same as 21 except for moving the peak to 0.61 is done properly, using 117, 
241, and 447 points for error tolerances 10 -3 , 10 -6 , and 10 -9 , respectively. 
DQUAD is clearly more efficient and robust, based on these test integrals, 
than any other ANQR tested except QSUBA. DQUAD is more robust than 
QSUBA, but is less efficient for integrals where QSUBS does not need to 
divide the interval. 

3.2 Parameter studies ( 1) 

Testing an ANQR on a "random" set of test integrals, while instructive 
and a good start, is insufficient for a library routine. The testing of an 
ANQR is incomplete without numerous parameter studies: 

.6 

f(x;a)dx 



s: 



with fixed e and varying a (Ref. 12), and with fixed a and varying e. The 
function f(x,a) should be increasingly difficult to integrate as a ap- 
proaches some limit, and a should be pushed close enough to that limit 
so that failure occurs. For error tolerance studies, the requested error 
should range from the approximate value of the integral to less than 
typical roundoff on the computer being used. 

Several of the first type of parameter study will now be discussed. For 
all of them, the error requested is 10 -6 . The first was suggested by de 
Boor 6 : 

2« 



X 



dx 



'o l + (2 a x) 2 

For large a, the integrand is highly peaked. The integrand is 2 a at x = 
0, falls to half that at x = 2~ a , and is 2 _a at the endpoints of the interval. 
There is no danger in missing the narrow spike, since it is exactly at the 
center of the interval, a normal sampling point for almost all ANQRs. This 
example demonstrates the power of adaptive ANQR, in that the number 
of sampling points increases only as \fa , approximately. The behavior 
of DQUAD is shown in Fig. 1. The number of sampling points is plotted 
against a. The meaning of the symbols used for plotting in all the figures 
is given below. 

• Successful integration; no error flag 
X Unsuccessful integration; error flag 
O Unsuccessful integration; no error flag 
® Successful integration; error flag 

In testing DQUAD, "successful" means that the true error is less than e. 
DQUAD failed for a > 38, but recognized its failure (failure was due to 
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Fig. 1 — Performance of DQUAD on function highly peaked at center of interval. 

filling of the function value stack). For comparison, CADRE fails starting 
at a = 31, and generally uses about 20 percent fewer sampling points than 
DQUAD. 

This integral should not pose any serious difficulty to a decent ANQR, 
since the unpleasant behavior occurred exactly at the center of the in- 
terval. Usually, though, the user should strive to break up integrals so 
that any unpleasant behavior happens at one of the endpoints of the 
interval. It is feasible for ANQRs to recognize such behavior at the end- 
points, but difficult if it occurs in the center of the interval. As an ex- 
ample, the previous integral, except with limits and 1, may be consid- 
ered. DQUAD's performance is shown in Fig. 2. For a ^ 7, DQUAD recog- 
nizes that the integrand approximates a step function, and makes a 
change of variable. This change of variable keeps the number of sampling 
points from growing significantly as a increases. 

Figure 3 shows a similar integral, except more steeply decaying away 
from the endpoint. 



i 
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Fig. 2 — DQUAD's performance on function highly peaked at end of interval. 
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Fig. 3 — Performance of DQUAD on steeply decaying function. 

Failure eventually occurs, as for the previous integral, when, after the 
change of variable, the new integrand is a sharply-peaked function with 
its peak away from an endpoint, and the peak is missed entirely. 

Testing of ANQRs on functions with isolated peaks within the range 
of integration takes more work; a must parameterize the narrowness of 
the peak, but the position of the peak is also important. A suitable test 
integral has two parameters. DQUAD was tested on 



£ 



2 « e -4«(*-0)2 dx 



For each a, 25 integrals were done, with 0's of 0.02(0.02.)0.50, for a sta- 
tistical evaluation. No failures occurred for a = 1, 2, 3, 4, 5, or 6; one oc- 
curred for a = 7, at = 0.04. For a = 8, 12 out of 25 failed. Figure 4 il- 
lustrates the results for = 0.40, a typical value. 

Another standard test integral, also used by de Boor, is 



£ 



x a dx 
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Fig. 4 — Performance of DQUAD on highly peaked function. 
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Fig. 5 — Performance of DQUAD on function with endpoint singularity. 

with the integrand set equal to zero at x = 0. Thus for a = the integrand 
is a step function, and for a < the integrand has an infinite disconti- 
nuity. For a < 0, the integrand is "unreasonable" by almost anyone's 
definition, but users sometimes give such integrals to ANQRs. Figure 5 
shows the performance of DQUAD. A change of variable was made au- 
tomatically by DQUAD for -0.97 < a < 1.75, approximately. DQUAD is 
designed to reject the change of variable if the estimated a is less than 
—0.99, and the change is not necessary to achieve the desired accuracy 
for a ^ 1.75. For a < —0.97, DQUAD fails with an error flag, using about 
900 function samples. Machine precision limits the efficiency of the 
change of variable for a less than about —%. For comparison, CADRE 
fails, with an error flag, for a less than — 7 /s> and gives an erroneous error 
flag for a near, but not at, a = 1. CADRE is substantially less efficient than 
DQUAD for this test integral. Other ANQRs, without special procedures 
for endpoint singularities, are much worse. 

A final type of test integral is one with an oscillatory integrand. Since 
almost all ANQRs sample the integrand at regularly spaced points, there 
is a danger of undersampling. As an example, consider 

• i 
[1 + cos (anx)] dx 



£ 
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Fig. 6 — Performance of DQUAD with oscillatory integrand. 

as suggested by de Boor. Figure 6 illustrates the performance of DQUAD 
on this test integral. The first failure is for a near 36, where the regular 
sampling causes the cosine to look like unity. The lower part of the figure 
shows the region near 36. The top part shows the number of sampling 
points used for a from V 4 to 20, with spacing of l / 4 . The trend line is ap- 
proximately N = 10a; CADRE's trend line is approximately N = 20a. 
Besides the general rising trend, there many dots significantly below the 
trend line. These occur because of resonance between the regular sam- 
pling of DQUAD and the regular oscillation of the integrand. If a is inte- 
gral, or very nearly so, coarse sampling may indicate that the integrand 
is simpler than it really is. For example, if a is an odd integer, the cosine 
is odd about x = 0.5, and the center third of the integral is integrated 
(correctly) based on insufficient sampling, since the trapezodial rule 
correctly integrates odd functions. 

This kind of resonance phenomenon is a difficult problem for an ANQR 
to surmount. The best way is probably to back off slightly from the no- 
tion of a fully automatic ANQR. Usually a prospective user of an ANQR 
knows if an integrand is oscillatory, and further, knows the approximate 
period of the oscillation. To avoid problems with possible undersampling, 
the user could divide the integral into several integrals, each with only 
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Fig. 7 — Parameter test with accuracy varied. 

a few periods of the oscillation. An easier way is for the user to require 
the ANQR to sample the integrand sufficiently finely to see all the os- 
cillations; this assumes the user knows the approximate period. DQUAD 
has provision for this, using the parameter HSAMPL, discussed in Section 
IV. If HSAMPL is taken to be \/a in the above example, all the dots fall 
on the main trend line, and no failures occur. 

3.3 Parameter studies (2) 

The second type of parameter testing, using a single test integral and 
varying the requested error, will be considered now. Some ANQRs will 
fail when the requested error is large compared to the value of the inte- 
gral, because of insufficiently cautious error estimation. All ANQRs fail 
if the requested error tolerance is too small, because of roundoff. A 
properly designed library ANQR will have some mechanism for dealing 
with too-small error tolerances. Finally, the graph of the number of 
sampling points versus the requested accuracy is of interest. 

Figure 7, for 



£ 



8e~ 8x dx 



is typical of DQUAD's performance on easy integrals. The graph of N 
against log 1/e is horizontal at large requested error, with no failures, since 
DQUAD does not believe any answer until it can be confident of the error 
estimate. The central portion of the graph is roughly linear. The number 
of sampling points needed is approximately proportional to c _1/2 °. 
DQUAD is effectively functioning as a 20th-order method. The small- 
est-error part of the graph is also horizontal; accuracy is limited by 
roundoff. 
Figure 8, for 

J x l ' 2 dx 

is similar except for the smallest requested errors. A change of variable 
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Fig. 8 — Parameter test with accuracy varied. 

is made only for e = 10 -4 through 10" 16 . For e = 10" 2 through 10 -15 , 
DQUAD is effectively a 15th-order method. The effective order is less than 
for the previous example because only three columns of the T table can 
be asymptotic. Roundoff begins to pollute the answer at e = 10~ 16 . 
Figure 9, for 



J[l + cos (airx)] dx 
o 



for two values of a, 1.95 and 17.95, is the last example. For a = 1.95, 
DQUAD is effectively a 20th-order method for e = 10" 3 through 10~ 17 . 
For a = 17.95, DQUAD does not begin to be a 20th-order method until 
e = 10" 6 . 

3.4 Parameter studies (3) 

The final type of parameter testing uses a single test integral and error 
tolerance, but varies the amount of noise in the function. For testing, 
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Fig. 9 — Parameter test with accuracy varied. 
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Fig. 10 — Parameter test with amount of noise varied. 



integrals of the types 



P [f(x) + 10 a r(x)]dx 
f /(*)[! + 10«r(x)] dx 



were done, with requested accuracy e = 10~ 6 . The function r(x) is the 
output of a pseudorandom number generator with values uniformly 
distributed on (-1,1). Values used for a were 1, 0, -1, . . . , -8. For/, the 
same four functions were used as in the previous section. For each 
function and each a, five different starting points of the random number 
generator were used. Sample results are shown in Fig. 10 and 11. The 
number of function values used is plotted vs. a, for each of the five tries. 
(Where two or more of the tries coincide, they are plotted side by side.) 
Points plotted with • are those for which DQUAD returned with an error 
estimate less than 10~ 6 ; points plotted with X indicate an error estimate 
greater than 10 -6 . 
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Fig. 11 — Parameter test with amount of noise varied. 
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A summary of the results follows. 

For every test with 10 a > 10 -6 , DQUAD recognized the presence of 
noise, estimated its magnitude, adjusted accuracy tolerances accordingly, 
and gave an answer and error estimate. The error estimates given were 
at most twice 10", and were generally less than 1.2 times 10". In only 6 
out of 280 integrals was the calculated answer farther from the noise-free 
value of the integral than the estimated error returned. 

For every test with 10" = 10" 6 , DQUAD returned an answer with an 
error estimate less than 10 -6 . In 3 out of 40 integrals, the presence of 
noise was recognized. In 2 out of 40, the calculated answer was between 
1 and 2 times 10 -6 away from the noise-free value; in the remaining 38, 
the answer was less than 10 -6 away. 

For every test with 10" < 10 -6 , DQUAD returned an answer with an 
error estimate of less than 10~ 6 , and the calculated answer was less than 
10~ 6 away from the noise-free value. No noise was recognized. 

IV. IMPLEMENTATION OF QUAD 

QUAD is a Fortran subroutine which attempts to evaluate 

>b f(x)dx 



I 



to within absolute accuracy t; the user supplies a, b, e, and a Fortran 
function subprogram to evaluate f(x). The discussion also applies to the 
double precision version, DQUAD, except as noted. 

QUAD is written in PFORT, 15 a portable subset of ANS Fortran. Tem- 
porary storage space is managed by a portable Fortran stack mecha- 
nism. 7 

The calling sequence is 

CALL QUAD (F,A,B,EPS,ANS,ERREST) 

F is the name of the user's subprogram, A and B are the limits, and EPS 
the requested accuracy. The estimated value of the integral and an es- 
timate of the accuracy of the estimate are given in ANS and ERREST. If 
ERREST is larger than EPS, QUAD invokes the PORT library's centralized 
error handling facility, 7 turning on the error state before returning 
control to the calling program. If the user has not taken prior action to 
recover from errors, an error message is printed and the run is termi- 
nated. The user who has taken prior action may test for the error state, 
and continue if desired. 

QUAD sets seven parameters to their default values and then calls 
RlQUAD. The default values should be adequate for almost all users, but 
the user wishing other values may call RlQUAD directly. RlQUAD also 
gives more information to the user about problems encountered during 
the integration. RlQUAD's calling sequence is 
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CALLRlQUAD(F,A,B,EPS,HSAMPL,NFCALL,LYSTAK,KMAXEX, 
KDIVID,JPRINT,NUMINT,ANS,ERREST,KWARN) 

The default parameter values used by QUAD are 



HSAMPL 


0.125 


NFCALL 


2000 


LYSTAK 


250 


KMAXEX 


6 


KDIVID 


4 


JPRINT 





NUMINT 


3 



The same parameters are used for DQUAD, except that KMAXEX = 8. 

HSAMPL measures how finely /(x) must be sampled. No error estimate 
is believed unless the trapezoidal rule step size is HSAMPL | b — a \ at 
smaller. Reducing HSAMPL improves the robustness of RlQUAD, but 
decreases its efficiency. Changing HSAMPL is useful for integrating os- 
cillatory functions, where there is some danger of aliasing. For Jo cos 
(airx ) dx, HSAMPL = 1/a is safe. The minimum number of sampling 
points of f(x) is roughly 2/hsampl. 

NFCALL is the upper limit on the number of sampling points of 
fix). 

LYSTAK is the length of the stack for storing values of the function 
at the sampling points. Space for the stack is allocated within RlQUAD 
using a portable storage -allocation package. 7 

KMAXEX is the maximum number of extrapolations allowed. Legal 
values are 4, 6, 8, 10, and 12. 

KDIVID is the minimum number of extrapolations required before 
dividing an interval. Legal values are 4, 6, 8, ... , KMAXEX. 

JPRINT determines how much printing will be done. With jprint < 
1, there is none. With JPRINT = 1, the endpoints, attempted error tol- 
erance, answer, and error estimate are printed for each interval at- 
tempted. With JPRINT > 1, the T tables are also printed. 

NUMINT is the number of equal intervals into which [a, b] is divided 
to start. If NUMINT > 1, signularities can be recognized at x = a and at 
x = b; if NUMINT = 1, a singularity can be recognized only at x = a. In- 
creasing NUMINT generally increases robustness while decreasing effi- 
ciency. 

KWARN is an integer warning flag, output from RlQUAD. It is zero if 
all went well. If KWARN is positive, it may have up to 6 digits, each with 
an independent meaning. Although ERREST may be greater than the 
requested EPS, it is still reliable. Each digit is zero unless a problem oc- 
curred; starting with the right-most digit, the problems are: 

1676 THE BELL SYSTEM TECHNICAL JOURNAL, NOVEMBER 1977 



(i) The error estimate is limited by noise or roundoff, but is above 
the requested error. 

(ii) The interval size is as small as is allowed. 

(Hi) As many intervals are stacked as is allowed. 

(iu) As many function values are stacked as is allowed. 

(u) As many function sampling points have been used as is al- 
lowed. 

(vi) The error estimate is larger than the requested error. 



4. 1 Machine-dependent constants 

All machine dependency is isolated into four machine-dependent 
constants, which are set by RlQUAD. No reprogramming is necessary to 
run QUAD or RlQUAD on another computer. The constants are set using 
the portable machine constants program RlMACH of Ref. 7. 

DLARGE is used as error estimate before any error estimates are con- 
sidered believable. Its value is set to slightly less than the largest float- 
ing-point number, 0.1 RlMACH(2). 

DSMALL is used as the default magnitude of r 2 (x). Its value is slightly 
larger than the smallest positive floating-point number, 10 
RlMACH(l). 

DROUND is used as the default magnitude of r\(x). Its value is set to 
50 times the largest relative rounding error, or 50 RlMACH(4). 

HSMALL is the smallest fraction of | b — a | used for a trapezoidal rule 
step size. Its value is the same as DROUND. 

REFERENCES 

1. F. Bauer, H. Rutishauser, and E. Stiefel, "New Aspects in Numerical Quadrature," 

Proc. of Symposia in Applied Mathematics, 15, pp. 199-218, American Mathe- 
matical Society, Providence, 1963. 

2. R. Bulirsch and J. Stoer, "Fehlerabschaetzungen and Extrapolation mit Rationalen 

Funktionen bei Verfahren vom Richardson-Typus," Num. Math., 6, 1964, pp. 
413-427. 

3. R. Bulirsch and J. Stoer, "Asymptotic Upper and Lower Bounds for Results of Ex- 

trapolation Methods," Num. Math., 8, 1966, pp. 93-104. 

4. W. Bunton, M. Diethelm, and K. Haigler, "Romberg Quadrature Schemes for Single 

and Multiple Integrals," Jet Propulsion Laboratory Report TM-324-221, 1969. 

5. P. J. Davis and P. Rabinowitz, Numerical Integration, Waltham, Mass.: Blaisdell, 

1967. 

6. C. de Boor, "CADRE: An Algorithm for Numerical Quadrature," in Mathematical 

Software, J. R. Rice (ed.), New York: Academic Press, 1971, pp. 417-449. 

7. P. A. Fox, A. D. Hall, and N. L. Schryer, "The PORT Mathematical Subroutine Li- 

brary," Bell Laboratories Computer Science Technical Report No. 47, September, 
1976. Accepted for publication in ACM Transactions on Mathematical Software. 
Inquiries about the PORT library may be addressed to Bell Laboratories Computing 
Information Service, 600 Mountain Avenue, Murray Hill, New Jersey 07974. 

8. D. K. Kahaner, "Comparison of Numerical Quadrature Formulas," in Mathematical 

Software, J. R. Rice (ed.), New York: Academic Press, 1971, pp. 229-259. 

9. R. E. Lynch, "Generalized Trapezoidal Formulas and Errors in Romberg Quadrature," 

Blanch anniversary volume, Aerospace Research Laboratories, Office of Aerospace 
Research, United States Air Force, pp. 215-229, 1967. 

AUTOMATIC NUMERICAL QUADRATURE 1677 



10. J. N. Lyness and B. W. Ninham, "Numerical Quadrature and Asymptotic Expansions," 

Math. Comp., 21 (1967), pp. 162-178. 

11. J. N. Lyness, "Algorithm 379, SQUANK," Comm. ACM, 13, 1970, pp. 260-263. 

12. J. N. Lyness, attribution by C. de Boor (Ref. 6). 

13. H. O'Hara and F. Smith, "The Evaluation of Definite Integrals by Interval Subdivi- 

sion," Nat. Bur. of Standards, Report N69- 11541, 1969. 

14. T. N. L. Patterson, "Algorithm 468, Algorithm for Automatic Numerical Integration 

over a Finite Interval," Comm. ACM, 16, 1973, pp. 694-699. 

15. B. G. Ryder, "The FORTRAN Verifier: User's Guide," Bell Telephone Laboratories 

Computing Science Technical Report No. 12, March 1973. 



1678 THE BELL SYSTEM TECHNICAL JOURNAL, NOVEMBER 1977 



